Motivation
Images can always contain a lot of information
Main character, objects, backgrounds, etc.

To help computer identify different elements in an image, we use segmentation to divide image into sub-regions.
Definition
Thresholding based segmentation directly divide images into regions based on intensity values and/or properties of these values.
Suppose that the intensity histogram in this plot corresponds to an image

To extract objects from background, one obvious way is to set a threshold that separates these modes.
Let f(x, y) being an image, and T being the threshold
If f(x, y) > T, these points are called object points
Otherwise, those points are called background points
In a segmented image, $g(x, y)= \begin{cases}1 & \text { if } f(x, y)>T \\ 0 & \text { if } f(x, y) \leq T\end{cases}$
There could also be some more difficult problems
In this case, the segmented image is given by $g(x, y)= \begin{cases}a & \text { if } f(x, y)>T_2 \\ b & \text { if } T_1 < f(x, y) \leq T_2 \\ c & \text { if } f(x, y) \leq T_1\end{cases}$
where a, b, and c are any three distinct intensity values.
Having noise in the image can influence the histogram of an image
Illumination and reflectance can also affect the histogram of an image and thus the quality of thresholding
Basic Global Thresholding:
Finding Optimum: Limitations and solution
Thresholding can be viewed as a statistical-decision theory problem.
Objective: Minimize the average error incurred in assigning pixels to two or more groups.
Elegant closed-form solution: Bayes decision rule
Solution is based on:
Issue: Estimating PDF is not a trivial task.
Alternative: Otsu's Method (Otsu [1979])
Why is it optimum?
Otsu's method maximizes the between-class variance, while variance is widely used in discriminant analysis.
Basic idea: Well-thresholded classes should be distinct with respect to the intensity values.
Therefore, such a threshold gives the best separation between classes.
Not only optimal, but also simple since entirely based on computations performed on this histogram.
The method can be summarized as follows:
import numpy as np
import os
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
import skimage as ski
from skimage.color import rgb2gray
from PIL import Image
import cv2
import random
file_name = os.path.join('griez.jpeg')
our_image = ski.io.imread(file_name)
img = mpimg.imread(file_name)
imgplot = plt.imshow(img)
gs_image = ski.color.rgb2gray(img)
plt.imshow(gs_image, cmap = 'gray')
<matplotlib.image.AxesImage at 0x7fd02ac13340>
n, bins, patches = plt.hist(gs_image.flatten(), bins=100)
var_list = []
for k in bins[1:]:
p1 = 0
p2 = 0
for i in range(len(bins)-1):
if bins[i] < k:
p1 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
else:
p2 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
mk = 0
mk2 = 0
for i in range(len(bins)-1):
if bins[i] < k:
mk += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
else:
mk2 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
var = p1*p2*(mk-mk2)**2
var_list.append(var)
ix = var_list.index(max(var_list))
print(bins[ix])
0.38
for i in range(len(gs_image[:,0])):
for j in range(len(gs_image[0,:])):
if gs_image[i][j] < bins[ix]:
gs_image[i][j] = 0
else:
gs_image[i][j] = 1
plt.imshow(gs_image, cmap = 'gray')
<matplotlib.image.AxesImage at 0x7fcfc86379d0>
Thresholding techniques discussed before can be extended to arbitrary number of classes.
In the case of $K$ classes, $C_1, C_2, ..., C_K$m the between class variance can be expressed as
$\sigma_B^2=\sum_{k=1}^K P_k(m_k-m_G)^2$
For example, if we have two threshold, the between class variance is
$\sigma_B^2=P_1(m_1-m_G)^2+P_2(m_2-m_G)^2+P_3(m_3-m_G)^2$
Iteratively performing this method over the entire image will give us two optimal thresholds
file_name = os.path.join('griez.jpeg')
our_image = ski.io.imread(file_name)
img = mpimg.imread(file_name)
gs_image = ski.color.rgb2gray(img)
mg = 0
for i in range(len(bins)):
if i < 100:
mg += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
var_array = np.zeros((len(bins[1:]), len(bins[1:])))
var_list2 = []
max_idx = []
for k in range(len(bins[1:])):
for k2 in range(len(bins[k+1:])):
p1 = 0
p2 = 0
p3 = 0
for i in range(len(bins)-1):
if i <= k:
p1 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
elif i <= k2:
p2 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
else:
p3 += n[i]/(np.shape(img)[0]*np.shape(img)[1])
mk1 = 0
mk2 = 0
mk3 = 0
for i in range(len(bins)-1):
if i <= k:
mk1 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
elif i <= k2:
mk2 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
else:
mk3 += bins[i]*n[i]/(np.shape(img)[0]*np.shape(img)[1])
var = p1 * (mk1-mg)**2 + p2 * (mk2-mg)**2 + p3 * (mk3-mg)**2
var_array[k,k2] = var
for i in range(len(var_array[0])):
max_idx.append(np.max(var_array[i]))
np.where(var_array == max(max_idx))
(array([38]), array([61]))
k1 = bins[38]
k2 = bins[61]
for i in range(len(gs_image[:,0])):
for j in range(len(gs_image[0,:])):
if gs_image[i][j] <= k1:
gs_image[i][j] = 0
elif k1 < gs_image[i][j] <= k2:
gs_image[i][j] = 0.5
else:
gs_image[i][j] = 1
plt.imshow(gs_image, cmap ='gray')
<matplotlib.image.AxesImage at 0x7fd0099f1700>
One of the simplest approach: Image partition
Example:
Methods:
Reference: [1]. Gonzalez, Rafael C., and Richard E. Woods. Digital Image Processing. Pearson, 2019. [2]. Mohammed, Zhana Fidakar, and Alan Anwer Abdulla. “Thresholding-based white blood cells segmentation from microscopic blood images.” UHD Journal of Science and Technology, vol. 4, no. 1, 13 Feb. 2020, pp. 9–17, https://doi.org/10.21928/uhdjst.v4n1y2020.pp9-17.